/*
 *	surface_normal
 *
 *	process the solid instrument geometry to find all the surface points and 
 *	compute the surface normal at each of these points
 *
 *	read the solid geometry from file geometry_instrument.dot produced by geom_3d_v4
 *
 *	there are two kinds of surface points:
 *		type A: which have one or more nearest neighbors (nn) open
 *		type B: which have one next nearest neighbors open - these are 
 *			concave ridge lines
 *		type C: which have one third nearest neighbors open - these are 
 *			concave corners 
 *
 *	usage: surface_normal geom_file surface_normal_file
 *		where geom_file is a geometry_instrument.dat file produced by geom_3d_v4
 */

/*
#define	N_X	400
#define N_Y	400
#define N_Z 	400
*/

#include <stdio.h>
#include <math.h>
#include <strings.h>
#include <string.h>
#include <ctype.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <omp.h>

#include "grid.h"

#define OPEN		0
#define TYPE_A		1
#define TYPE_B		2
#define TYPE_C		3
#define INTERIOR	4

#define N_S	2000000

char point_type[N_X_MAX][N_Y_MAX][N_Z_MAX];

int n_surface;

int i_surface[N_S];
int j_surface[N_S];
int k_surface[N_S];
double i_hat[N_S];
double j_hat[N_S];
double k_hat[N_S];

int i_min,i_max,j_min,j_max,k_min,k_max;

int n_A,n_B,n_C;

main(int argc,char *argv[])
{
	if(argc != 2) {
		fprintf(stderr,"usage: surface_normal geom_file\n");
		(void)exit(0);
	}
	init_grid();
	readin(argv[1]);
fprintf(stderr,"back from readin\n");
	find_surface_points_A();
fprintf(stderr,"back from find surface points\n");
	find_surface_points_B();
	save_surface_normals();
}

init_grid()
{
	int i,j,k;

	for(i = 0; i < N_X_MAX; i++) {
		for(j = 0; j < N_Y_MAX; j++) {
			for(k = 0; k < N_Z_MAX; k++) {
				point_type[i][j][k] = OPEN;
			}
		}
	}

	return;
}

save_surface_normals()
{
	int i,j,k;
	FILE *fp_out;
	double norm;

// 	put info on suface points and their surface normals into file surface_point_info.dat
//	store as i,j,k coordinates followed by i_hat,j_hat,k_hat
//	note that i_hat, etc., are doubles (not ints)

	fp_out = fopen("surface_point_info.dat","w");

	for(i = 0; i < n_surface; i++) {
		fprintf(fp_out,"%d\t%d\t%d\t%g\t%g\t%g\n",i_surface[i],j_surface[i],
			k_surface[i],i_hat[i],j_hat[i],k_hat[i]);
// printf("%d %g %g %g\n",i,i_hat[i],j_hat[i],k_hat[i]);
	}
	fclose(fp_out);
	
	return;
}

find_surface_points_A()
{
	int i,j,k;
	double  x_sum,y_sum,z_sum;
	int n_open;
	double norm;

	int n_nearest_neighbors; // number of occupied nn

//	look for type A points
	n_A = 0;
	n_surface = 0;

fprintf(stderr,"starting find surface points\n");
fprintf(stderr,"%d\t%d\t%d\t%d\t%d\t%d\n",i_min,i_max,j_min,j_max,k_min,k_max);

	for(i = i_min; i <= i_max; i++) {
		for(j = j_min; j <= j_max; j++) {
			for(k = k_min; k <= k_max; k++) {
			   if(point_type[i][j][k] == INTERIOR) {
				n_nearest_neighbors = 0; 
				x_sum = y_sum = z_sum = 0;
				n_open = 0;
				if(point_type[i+1][j][k] == OPEN) { 
					x_sum += 1;
					++n_open;
				}
				if(point_type[i-1][j][k] == OPEN) {
					x_sum -= 1;
					++n_open;
				}
				if(point_type[i][j+1][k] == OPEN) {
					y_sum += 1;
					++n_open;
				}
				if(point_type[i][j-1][k] == OPEN) {
					y_sum -= 1;
					++n_open;
				}
				if(point_type[i][j][k+1] == OPEN) {
					z_sum += 1;
					++n_open;
				}
				if(point_type[i][j][k-1] == OPEN) {
					z_sum -= 1;
					++n_open;
				}
		
				if(n_open > 0) {
					point_type[i][j][k] = TYPE_A;
					i_surface[n_surface] = i;
					j_surface[n_surface] = j;
					k_surface[n_surface] = k;
					i_hat[n_surface] = x_sum / n_open;
					j_hat[n_surface] = y_sum / n_open;
					k_hat[n_surface] = z_sum / n_open;
					norm = i_hat[n_surface]*i_hat[n_surface] 
						+ j_hat[n_surface]*j_hat[n_surface] 
						+ k_hat[n_surface]*k_hat[n_surface];
					norm = sqrt(norm);
					i_hat[n_surface] /= norm;
					j_hat[n_surface] /= norm;
					k_hat[n_surface] /= norm;
// printf("%d %d %d %g %g %g\n",i,j,k,i_hat[n_surface],j_hat[n_surface],k_hat[n_surface]);
					++n_surface;
					++n_A;
				}	
			   }
			}
		}
	}
	fprintf(stderr,"number of type A surface points = %d\n",n_A);
	fprintf(stderr,"number of surface points = %d\n",n_surface);

	return;
}

int readin(char *s)
{
	int i,j,k,n;
	FILE *fp_in;

	if((fp_in = fopen(s,"r")) == NULL) {
		fprintf(stderr,"can't open %s\n",s);
		(void)exit(0);
	}

	n = 0;
	if(fscanf(fp_in,"%d %d %d",&i,&j,&k) == EOF) {
		fprintf(stderr,"no data in %s\n",s);
		(void)exit(0);
	}
	i_min = i_max = i;
	j_min = j_max = j;
	k_min = k_max = k;
	n = 1;

	while(1) {
		if(fscanf(fp_in,"%d %d %d",&i,&j,&k) == EOF) break;
		if(i_min > i) i_min = i;
		if(i_max < i) i_max = i;

		if(j_min > j) j_min = j;
		if(j_max < j) j_max = j;

		if(k_min > k) k_min = k;
		if(k_max < k) k_max = k;
		++n;
	}
	fclose(fp_in);

fprintf(stderr,"n points = %d\n",n);

	fp_in = fopen(s,"r");

	while(1) {
		if(fscanf(fp_in,"%d %d %d",&i,&j,&k) == EOF) break;
		point_type[i][j][k] = INTERIOR;
	}
	fclose(fp_in);

	return;
}
	
find_surface_points_B()
{
	int i,j,k;
	double  x_sum,y_sum,z_sum;
	int n_open;
	double norm;

	int n_next_nearest_neighbors; // number of occupied nn

//	look for type B points
	n_B = 0;

fprintf(stderr,"starting find surface points B\n");

	for(i = i_min; i <= i_max; i++) {
		for(j = j_min; j <= j_max; j++) {
			for(k = k_min; k <= k_max; k++) {
			   if(point_type[i][j][k] == INTERIOR) {
				n_next_nearest_neighbors = 0; 
				x_sum = y_sum = z_sum = 0;
				n_open = 0;
				if(point_type[i+1][j+1][k] == OPEN) { 
					x_sum += 1;
					y_sum += 1;
					++n_open;
				}
				if(point_type[i+1][j-1][k] == OPEN) { 
					x_sum += 1;
					y_sum -= 1;
					++n_open;
				}
				if(point_type[i+1][j][k+1] == OPEN) { 
					x_sum += 1;
					z_sum += 1;
					++n_open;
				}
				if(point_type[i+1][j][k-1] == OPEN) { 
					x_sum += 1;
					z_sum -= 1;
					++n_open;
				}
				if(point_type[i-1][j+1][k] == OPEN) {
					x_sum -= 1;
					y_sum += 1;
					++n_open;
				}
				if(point_type[i-1][j][k+1] == OPEN) {
					x_sum -= 1;
					z_sum += 1;
					++n_open;
				}
				if(point_type[i-1][j-1][k] == OPEN) {
					x_sum -= 1;
					y_sum -= 1;
					++n_open;
				}
				if(point_type[i-1][j][k-1] == OPEN) {
					x_sum -= 1;
					z_sum -= 1;
					++n_open;
				}
				if(point_type[i][j+1][k+1] == OPEN) {
					y_sum += 1;
					z_sum += 1;
					++n_open;
				}
				if(point_type[i][j+1][k-1] == OPEN) {
					y_sum += 1;
					z_sum -= 1;
					++n_open;
				}
				if(point_type[i][j-1][k+1] == OPEN) {
					y_sum -= 1;
					z_sum += 1;
					++n_open;
				}
				if(point_type[i][j-1][k-1] == OPEN) {
					y_sum -= 1;
					z_sum -= 1;
					++n_open;
				}
		
				if(n_open > 0) {
					point_type[i][j][k] = TYPE_B;
					i_surface[n_surface] = i;
					j_surface[n_surface] = j;
					k_surface[n_surface] = k;
					i_hat[n_surface] = x_sum / n_open;
					j_hat[n_surface] = y_sum / n_open;
					k_hat[n_surface] = z_sum / n_open;
					norm = i_hat[n_surface]*i_hat[n_surface] 
						+ j_hat[n_surface]*j_hat[n_surface] 
						+ k_hat[n_surface]*k_hat[n_surface];
					norm = sqrt(norm);
					i_hat[n_surface] /= norm;
					j_hat[n_surface] /= norm;
					k_hat[n_surface] /= norm;
// printf("%d %d %d %g %g %g\n",i,j,k,i_hat[n_surface],j_hat[n_surface],k_hat[n_surface]);
					++n_surface;
					++n_B;
				}	
			   }
			}
		}
	}
	fprintf(stderr,"number of type B surface points = %d\n",n_B);
	fprintf(stderr,"total number of surface points = %d\n",n_surface);

	return;
}
